#install.packages("/Users/yazhuodeng/Downloads/nlmeU_0.70-3.tar", repos = NULL, type = "source")

library(nlmeU)
## 
## Attaching package: 'nlmeU'
## The following object is masked from 'package:stats':
## 
##     sigma
#data(package = "nlmeU")

library(lme4)   # using lmer()
## Loading required package: Matrix
library(lattice)
library(nlme)
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:lme4':
## 
##     lmList

Progressive Resistance Training (PRT) Trial

subjects’ characteristics: prt.subjects

fiber measurement: prt.fiber

merge subjects and fiber measurement: prt

str(prt.subjects)
## 'data.frame':    63 obs. of  5 variables:
##  $ id   : Factor w/ 63 levels "5","10","15",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ prt.f: Factor w/ 2 levels "High","Low": 2 1 1 1 1 2 2 2 2 1 ...
##  $ age.f: Factor w/ 2 levels "Young","Old": 1 1 2 2 2 1 1 2 1 1 ...
##  $ sex.f: Factor w/ 2 levels "Female","Male": 1 1 1 1 1 1 2 1 2 2 ...
##  $ bmi  : num  25.4 22.6 25.4 25.8 27.7 ...
str(prt.fiber)
## 'data.frame':    2471 obs. of  5 variables:
##  $ id     : Factor w/ 63 levels "5","10","15",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ iso.fo : num  0.265 0.518 0.491 0.718 0.16 0.41 0.371 0.792 0.649 0.604 ...
##  $ spec.fo: num  83.5 132.8 161.1 158.8 117.9 ...
##  $ occ.f  : Factor w/ 2 levels "Pre","Pos": 1 1 1 1 1 1 1 1 1 1 ...
##  $ fiber.f: Factor w/ 2 levels "Type 1","Type 2": 1 1 2 1 2 1 1 1 2 1 ...
str(prt)
## 'data.frame':    2471 obs. of  9 variables:
##  $ id     : Factor w/ 63 levels "5","10","15",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ prt.f  : Factor w/ 2 levels "High","Low": 2 2 2 2 2 2 2 2 2 2 ...
##  $ age.f  : Factor w/ 2 levels "Young","Old": 1 1 1 1 1 1 1 1 1 1 ...
##  $ sex.f  : Factor w/ 2 levels "Female","Male": 1 1 1 1 1 1 1 1 1 1 ...
##  $ bmi    : num  25.4 25.4 25.4 25.4 25.4 ...
##  $ iso.fo : num  0.265 0.518 0.491 0.718 0.16 0.41 0.371 0.792 0.649 0.604 ...
##  $ spec.fo: num  83.5 132.8 161.1 158.8 117.9 ...
##  $ occ.f  : Factor w/ 2 levels "Pre","Pos": 1 1 1 1 1 1 1 1 1 1 ...
##  $ fiber.f: Factor w/ 2 levels "Type 1","Type 2": 1 1 2 1 2 1 1 1 2 1 ...
#summary stats for subjects' characteristics
with(prt.subjects, tapply(bmi,prt.f, summary))
## $High
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.37   22.86   24.77   25.11   28.19   31.00 
## 
## $Low
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.95   23.05   24.84   24.70   26.27   32.29
by(subset(prt.subjects, select = -id), prt.subjects$prt.f, summary)
## prt.subjects$prt.f: High
##   prt.f      age.f       sex.f         bmi       
##  High:31   Young:15   Female:17   Min.   :18.37  
##  Low : 0   Old  :16   Male  :14   1st Qu.:22.86  
##                                   Median :24.77  
##                                   Mean   :25.11  
##                                   3rd Qu.:28.19  
##                                   Max.   :31.00  
## -------------------------------------------------------- 
## prt.subjects$prt.f: Low
##   prt.f      age.f       sex.f         bmi       
##  High: 0   Young:15   Female:17   Min.   :18.95  
##  Low :32   Old  :17   Male  :15   1st Qu.:23.05  
##                                   Median :24.84  
##                                   Mean   :24.70  
##                                   3rd Qu.:26.27  
##                                   Max.   :32.29
# number of fibers per type and occasion for the subject "5" and "335"
fibL<-
  with(prt,
       tapply(spec.fo,
              list(id = id, fiberF = fiber.f, occF = occ.f),
              length))
dimnms<- dimnames(fibL)
names(dimnms)
## [1] "id"     "fiberF" "occF"
fibL["5", , ]
##         occF
## fiberF   Pre Pos
##   Type 1  12  18
##   Type 2   7   4
fibL["335", , ]
##         occF
## fiberF   Pre Pos
##   Type 1  NA   8
##   Type 2  14  11
# mean value of spec.fo by fiber type and occasion for subject 5
fibM<-
  with(prt,
       tapply(spec.fo,
              list(id = id, fiberF = fiber.f, occF = occ.f),
              mean))
fibM["5", , ]
##         occF
## fiberF        Pre      Pos
##   Type 1 132.5917 129.9556
##   Type 2 145.7429 147.9500
#end at page 51

#  (a) Preprocessing of the data (melting)





#  (b) Aggregating data (casting)

Study of Instructional Improvement (SII)

variable information page 25

str(SIIdata)
## 'data.frame':    1190 obs. of  12 variables:
##  $ sex     : Factor w/ 2 levels "M","F": 2 1 2 1 1 2 1 1 2 1 ...
##  $ minority: Factor w/ 2 levels "Mnrt=No","Mnrt=Yes": 2 2 2 2 2 2 2 2 2 2 ...
##  $ mathkind: int  448 460 511 449 425 450 452 443 422 480 ...
##  $ mathgain: int  32 109 56 83 53 65 51 66 88 -7 ...
##  $ ses     : num  0.46 -0.27 -0.03 -0.38 -0.03 0.76 -0.03 0.2 0.64 0.13 ...
##  $ yearstea: num  1 1 1 2 2 2 2 2 2 2 ...
##  $ mathknow: num  NA NA NA -0.11 -0.11 -0.11 -0.11 -0.11 -0.11 -0.11 ...
##  $ housepov: num  0.082 0.082 0.082 0.082 0.082 0.082 0.082 0.082 0.082 0.082 ...
##  $ mathprep: num  2 2 2 3.25 3.25 3.25 3.25 3.25 3.25 3.25 ...
##  $ classid : Factor w/ 312 levels "1","2","3","4",..: 160 160 160 217 217 217 217 217 217 217 ...
##  $ schoolid: Factor w/ 107 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ childid : Factor w/ 1190 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
#examine the data hierarchy 
dtId<-subset(SIIdata, select = c(schoolid, classid, childid))
any(duplicated(dtId))
## [1] FALSE
require(nlme)
names(gsummary(dtId, form = ~childid, inv = T))
## [1] "schoolid" "classid"  "childid"
names(gsummary(dtId, form = ~classid, inv = T))
## [1] "schoolid" "classid"
names(gsummary(dtId, form = ~schoolid, inv = T))
## [1] "schoolid"
#school level variables
(nms1<-names(gsummary(SIIdata, form = ~schoolid, inv = T)))
## [1] "housepov" "schoolid"
#Class level variables
nms2a<-names(gsummary(SIIdata, form = ~classid, inv = T))
idx1<-match(nms1,nms2a)
(nms2<-nms2a[-idx1])
## [1] "yearstea" "mathknow" "mathprep" "classid"
#pupil level variables
nms3a<-names(gsummary(SIIdata, form = ~childid, inv = T))
idx12<-match(c(nms1,nms2),nms3a)
nms3a[-idx12]
## [1] "sex"      "minority" "mathkind" "mathgain" "ses"      "childid"
#The number of missing values for variables included in the SIIdata data frame  
sapply(SIIdata, FUN = function(x) any(is.na(x)))
##      sex minority mathkind mathgain      ses yearstea mathknow housepov 
##    FALSE    FALSE    FALSE    FALSE    FALSE    FALSE     TRUE    FALSE 
## mathprep  classid schoolid  childid 
##    FALSE    FALSE    FALSE    FALSE
sum(as.numeric(is.na(SIIdata$mathknow)))
## [1] 109
range(SIIdata$mathknow, na.rm = TRUE)
## [1] -2.50  2.61
#Extracting the information about the number of pupils per school
(schlN <- xtabs(~schoolid, SIIdata)) # Number of pupils per school
## schoolid
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18 
##  11  10  14   6   6  12  14  16   6  18  31  27   9  15  13   6   9   4 
##  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36 
##  11  12  17   4   5   8   7  15  21  10   8   3  22   9  22   7  11   8 
##  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54 
##  11   9  12   7   9  18   4  17   5  13  10   6   8  12   4   9   8  10 
##  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72 
##  13  10   8   5   4  13  11  13   3  10   9   8  12  16   5  19  27  11 
##  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90 
##  10   3  15  27  24  15  12   7   9  19  12  14  20   7  10   6   4   3 
##  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 
##  16   9  21  15   5   8   6   2  19  13  16  11   8   6  10   2  10
range(schlN)
## [1]  2 31
xtabs(~schlN)  # distribution of the number of pupils over schools
## schlN
##  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 24 27 31 
##  2  4  6  5  8  5  9  9 10  7  7  6  3  5  4  2  2  3  1  2  2  1  3  1
# Computation of the mean value of pupils’ math scores for each school
attach(SIIdata)
mthgM <- by(cbind(mathgain, mathkind), schoolid, colMeans)
head(mthgM, 10)
## $`1`
##  mathgain  mathkind 
##  59.63636 458.36364 
## 
## $`2`
## mathgain mathkind 
##     65.0    487.9 
## 
## $`3`
##  mathgain  mathkind 
##  88.85714 469.14286 
## 
## $`4`
##  mathgain  mathkind 
##  35.16667 462.66667 
## 
## $`5`
##  mathgain  mathkind 
##  60.16667 445.33333 
## 
## $`6`
##  mathgain  mathkind 
##  78.66667 457.08333 
## 
## $`7`
##  mathgain  mathkind 
##  57.35714 442.07143 
## 
## $`8`
## mathgain mathkind 
##  47.2500 498.6875 
## 
## $`9`
##  mathgain  mathkind 
##  34.83333 475.33333 
## 
## $`10`
##  mathgain  mathkind 
##  74.94444 441.33333
detach(SIIdata)

#  Constructing a data frame with summary data for schools
# (a) Creating a data frame with the number of classes and children for each school
library(reshape)
## 
## Attaching package: 'reshape'
## The following object is masked from 'package:Matrix':
## 
##     expand
idvars <- c("schoolid")
mvars <- c("classid", "childid")
dtm1 <- melt(SIIdata, id.vars = idvars, measure.vars = mvars)  
names(cst1 <-cast(dtm1,
                  fun.aggregate = function(el) length(unique(el))))
## [1] "schoolid" "classid"  "childid"
names(cst1) <- c("schoolid", "clssn", "schlN")

# (b) Creating a data frame with the school-specific means of selected variables
mvars <- c("mathgain", "mathkind", "housepov")
dtm2 <- melt(SIIdata, id.vars = idvars, measure.vars = mvars)
names(cst2 <- cast(dtm2, fun.aggregate = mean))
## [1] "schoolid" "mathgain" "mathkind" "housepov"
names(cst2) <- c("schoolid", "mthgMn", "mthkMn", "housepov")

#(c) Merging the data frames created in parts (a) and (b) above
schlDt <- merge(cst1, cst2, sort = FALSE)
head(schlDt,15)
##    schoolid clssn schlN   mthgMn   mthkMn housepov
## 1         1     2    11 59.63636 458.3636    0.082
## 2         2     3    10 65.00000 487.9000    0.082
## 3         3     4    14 88.85714 469.1429    0.086
## 4         4     2     6 35.16667 462.6667    0.365
## 5         5     1     6 60.16667 445.3333    0.511
## 6         6     3    12 78.66667 457.0833    0.044
## 7         7     4    14 57.35714 442.0714    0.148
## 8         8     3    16 47.25000 498.6875    0.085
## 9         9     3     6 34.83333 475.3333    0.537
## 10       10     4    18 74.94444 441.3333    0.346
## 11       11     9    31 50.45161 499.8065    0.101
## 12       12     5    27 52.40741 462.5185    0.106
## 13       13     2     9 47.66667 479.4444    0.055
## 14       14     3    15 44.33333 468.9333    0.157
## 15       15     5    13 75.00000 445.9231    0.121
rm(cst1, cst2)

# Summary statistics for the school-specific mean values of housepov
summary(schlDt$housepov)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0120  0.0855  0.1480  0.1941  0.2645  0.5640
xyplot(mthgMn ~ housepov, 
       schlDt, type = c("p", "smooth"), grid = TRUE) # Fig.3.8a

xyplot(mthgMn ~ mthkMn, 
       schlDt, type = c("p", "smooth"), grid = TRUE) # Fig.3.8b

# Extracting the information about the number of pupils per class

(clssN <- xtabs(~ classid, SIIdata))
## classid
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18 
##   5   3   3   6   1   5   1   4   3   2   4   5   9   4   1   6   3   1 
##  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36 
##   2   3   4   8   1   1   4  10   6   5   3   4   4   4   4   1   4   1 
##  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54 
##   1   7   2   3   5  10   3   6   5   5   1   4   6   4   6   2   5   4 
##  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72 
##   5   5   4   1   4   3   2   2   2   5   1   6   2   4   3   1   2   3 
##  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90 
##   3   1   6   4   1   8   5   6   8   5   7   1   7   5   1   2   3   2 
##  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108 
##   4   1   1   6   2   8   2   5   1   1   2   3   3   6   4   2   4   3 
## 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 
##   4   4   3   2   3   4   4   2   8   1   1   8   6   6   8   1   4   6 
## 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 
##   3   4   3   4   2   3   3   4   7   6   2   5   4   2   5   4   6   3 
## 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 
##   5   3   5   4   4   1   2   6   6   2   3   6   4   8   4   3   4   2 
## 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 
##   1   5   6   4   5   7   6   2   1   4   2   7   3   4   3   2   2   6 
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 
##   1   5   3   5   1   6   1   4   9   8   1   2   2   2   4   3   2   2 
## 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 
##   1   6   5   3   2   2   9   5   6   4   2   3   4   4   3   1   2   2 
## 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 
##   8   2   4   2   2   4   4   5   3   2   2   3   5   2   8   2   1   7 
## 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 
##   6   2   3   1   2   3   5   7   4   1   3   4   5   6   1   2   3   6 
## 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 
##   9   2   6   1   4   3   8   3   1   5   3   7   4   5   3   5   2   3 
## 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 
##   3   1   5   1   3   7   7   8   7   5   4   2   2   3   3   6   5   4 
## 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 
##   5   3   4   4   4   3   1   7   5   5   6   5   2   7   4   4   4   4 
## 307 308 309 310 311 312 
##   4   3   3   3   2   4
sum(clssN) # Total number of pupils
## [1] 1190
range(clssN)
## [1]  1 10
(clssCnt <- xtabs(~clssN)) # Distribution of no. of pupils/classes
## clssN
##  1  2  3  4  5  6  7  8  9 10 
## 42 53 53 61 39 31 14 13  4  2
sum(clssCnt) # Total number of classes
## [1] 312

Model specification

Model M18.6 is defined as follows: \[ MATHGAIN_{sci}=\beta_0 + \beta_1*SES_{sci}+\beta_2*MINORITY_{sci} +\beta_{1,2}*SES_{sci}*MINORITY_{sci} \] \[ +\beta_{3,p1}* p1(MATHKIND_{sci})+\beta_{3,p2}* p2(MATHKIND_{sci})+\beta_{3,p3}* p3(MATHKIND_{sci}) \] \[ +b_{0s}+b{0sc}+ \epsilon_{sci} \equiv \mu_{sci}+b_{0s}+b{0sc}+ \epsilon_{sci} \] the random-effects structure for a two-level LMM with nested effects. In particular, the nesting of grouping factors, i.e., classid within schoolid, is explicitly expressed using the crossing operator : (see Sect. 5.2.1) in the Z-term (1|schoolid:classid), included in the model formula along with the (1|schoolid) term.

library(lme4)

fm18.6mer <- lmer(mathgain ~ ses + minority + poly(mathkind, 3) + 
                    ses:minority + (1|schoolid) + 
                    (1|schoolid:classid), # Syntax #1 (general)
                  data = SIIdata, REML = FALSE) 
summ <- summary(fm18.6mer)
print(summ, corr = FALSE)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: mathgain ~ ses + minority + poly(mathkind, 3) + ses:minority +  
##     (1 | schoolid) + (1 | schoolid:classid)
##    Data: SIIdata
## 
##      AIC      BIC   logLik deviance df.resid 
##  11347.8  11398.6  -5663.9  11327.8     1180 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.8391 -0.6355 -0.0174  0.5937  4.0555 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  schoolid:classid (Intercept)  85.98    9.273  
##  schoolid         (Intercept)  66.77    8.171  
##  Residual                     690.89   26.285  
## Number of obs: 1190, groups:  schoolid:classid, 312; schoolid, 107
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)            61.354      2.069  29.649
## ses                     8.845      1.942   4.554
## minorityMnrt=Yes       -6.859      2.280  -3.008
## poly(mathkind, 3)1   -660.458     30.989 -21.312
## poly(mathkind, 3)2    124.455     28.240   4.407
## poly(mathkind, 3)3   -178.336     27.975  -6.375
## ses:minorityMnrt=Yes   -5.820      2.438  -2.387

In the data frame SIIdata, the levels of the factor classid have been coded as explicitly nested within the levels of schoolid (Sect. 2.4.3). This approach is actually recommended in the case of representing factors with nested levels. Hence, it is possible to fit model M18.6 using a simpler syntax, namely, (1|schoolid) + (1|classid) for the random-effects structure (see syntax (2a) in Table 15.2). In particular, in the second syntax shown in Panel R18.17, the Z-term for the factor classid, (1|classid), does not use the crossing operator : and, therefore, does not explicitly indicate the nesting.

update(fm18.6mer,mathgain ~ ses + minority + 
         poly(mathkind, 3) + ses:minority +
         (1|schoolid) + (1|classid))
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: mathgain ~ ses + minority + poly(mathkind, 3) + (1 | schoolid) +  
##     (1 | classid) + ses:minority
##    Data: SIIdata
##      AIC      BIC   logLik deviance df.resid 
## 11347.82 11398.64 -5663.91 11327.82     1180 
## Random effects:
##  Groups   Name        Std.Dev.
##  classid  (Intercept)  9.273  
##  schoolid (Intercept)  8.171  
##  Residual             26.285  
## Number of obs: 1190, groups:  classid, 312; schoolid, 107
## Fixed Effects:
##          (Intercept)                   ses      minorityMnrt=Yes  
##               61.354                 8.845                -6.859  
##   poly(mathkind, 3)1    poly(mathkind, 3)2    poly(mathkind, 3)3  
##             -660.458               124.455              -178.336  
## ses:minorityMnrt=Yes  
##               -5.820

Finally, the third form of the lmer()-function call, shown in Panel R18.17, uses the nonessential operator / (see Table 5.3) in the Z-term (1|schoolid/classid) to abbreviate the specification of the random-effects part of the lmer() model formula.

update(fm18.6mer,mathgain ~ ses + minority + 
         poly(mathkind, 3) + ses:minority +
         (1|schoolid/classid))
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: 
## mathgain ~ ses + minority + poly(mathkind, 3) + (1 | schoolid/classid) +  
##     ses:minority
##    Data: SIIdata
##      AIC      BIC   logLik deviance df.resid 
## 11347.82 11398.64 -5663.91 11327.82     1180 
## Random effects:
##  Groups           Name        Std.Dev.
##  classid:schoolid (Intercept)  9.273  
##  schoolid         (Intercept)  8.171  
##  Residual                     26.285  
## Number of obs: 1190, groups:  classid:schoolid, 312; schoolid, 107
## Fixed Effects:
##          (Intercept)                   ses      minorityMnrt=Yes  
##               61.354                 8.845                -6.859  
##   poly(mathkind, 3)1    poly(mathkind, 3)2    poly(mathkind, 3)3  
##             -660.458               124.455              -178.336  
## ses:minorityMnrt=Yes  
##               -5.820
#Extracting information about the estimated fixed- and random- effects structure of model M18.6 from the mer-class model-fit object. 
anova(fm18.6mer) # Approximate F-test statistics
## Analysis of Variance Table
##                   Df Sum Sq Mean Sq  F value
## ses                1    482     482   0.6970
## minority           1      8       8   0.0119
## poly(mathkind, 3)  3 368145  122715 177.6181
## ses:minority       1   3938    3938   5.6992
logLik(fm18.6mer)  #ML value
## 'log Lik.' -5663.91 (df=10)
unlist(VarCorr(fm18.6mer))
## schoolid:classid         schoolid 
##         85.98146         66.76779
as.data.frame(VarCorr(fm18.6mer))[3,"sdcor"]   #sigma
## [1] 26.28482
#(a) Normal Q-Q plot of the raw class-level residuals
rsd6 <- resid(fm18.6mer)
qqnorm(rsd6)

#(b) Normal Q-Q plot of predicted random effects
rnf6mer <- ranef(fm18.6mer)     #Random effects
rnf6qn <- plot(rnf6mer, grid = TRUE)# Q-Q plot for random effects
update(rnf6qn[["schoolid:classid"]], # For classid (see Fig.18.9a)
       ylab = c("Random effects for classes")) 

update(rnf6qn[["schoolid"]], # For schoolid (see Fig.18.9b)
       ylab = c("Random effects for schools"))

Flemish Community Attainment-Targets (FCAT) Study

#page 35
str(fcat)
## 'data.frame':    4851 obs. of  3 variables:
##  $ target: Factor w/ 9 levels "T1(4)","T2(6)",..: 1 2 3 4 5 6 7 8 9 1 ...
##  $ id    : Factor w/ 539 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ scorec: int  4 6 4 1 7 6 6 5 5 3 ...
#Investigation of the data grouping structure
tab1 <- xtabs(~ id + target, data = fcat)   # id by target table

all(tab1 > 0) # All counts > 0? 
## [1] TRUE
range(tab1) # Range of counts
## [1] 1 1
#(a) Summarizing scores for each child and attainment target
scM <- with(fcat, tapply(scorec, list(id, target), mean))

#(b) Histograms of scores for different attainment targets

histogram(~scorec | target, data = fcat, breaks = NULL) # Fig.3.11